library(reshape2)
library(ggplot2)

######Load data 
jan15 <- readRDS("jan15.rds") ## January 15
mar15 <- readRDS("mar15.rds") ## March 15
nov20 <- readRDS("nov20.rds") ## November 2020
feb21 <- readRDS("feb21.rds") ## February 2021
mar21 <- readRDS("mar21.rds") ## March 2021
jun21 <- readRDS("jun21.rds") ## June 2021
jun22 <- readRDS("jun22.rds") ## June 2022

#####Figure 1
###############Merge direct estimates across waves
p.jun21<-apply(jun21[,c("ziuganov", "mironov", "zhirinovskii", "putin", "mandela",
                        "merkel", "nazarbaev", "stalin", "brezhnev", 
                        "gorbachev", "yeltsin", "mikhalkov", "grudinin",
                        "sobchak", "kudrin")], 2, mean, na.rm=T)
p.mar15<-apply(mar15[,c("ziuganov", "mironov", "zhirinovskii", "putin", 
                      "mandela", "merkel", "lukashenko", "stalin", 
                      "brezhnev", "yeltsin", "castro")], 2, mean, na.rm=T)
p.jan15<-apply(jan15[,c("ziuganov", "mironov", "zhirinovskii", "putin",
                      "stalin", "brezhnev", "yeltsin")], 2, 
               mean, na.rm=T)
p.nov20<-mean(nov20[,c("putin")], na.rm=T)
p.mar21<-apply(mar21[,c("ziuganov", "mironov", "zhirinovskii", "putin", "mandela",
                        "merkel", "nazarbaev", "stalin", "brezhnev", 
                        "yeltsin", "mikhalkov", "grudinin",
                        "sobchak", "castro","navalny")], 2, mean, na.rm=T)
p.feb21<-apply(feb21[,c("ziuganov", "mironov", "zhirinovskii", "putin",
                        "stalin", "brezhnev", 
                        "yeltsin", "mikhalkov", "grudinin",
                        "sobchak", "navalny")], 2, mean, na.rm=T)
p.jun22<-apply(jun22[,c("putin",
                        "stalin", "brezhnev", "yeltsin", "castro")], 2, 
               mean, na.rm=T)

dir.pop<-data.frame(matrix(ncol=8,nrow=18))
colnames(dir.pop)<-c("Figure","2015-01-01","2015-03-01",
                     "2020-11-01", "2021-02-01", "2021-03-01",
                     "2021-06-01","2022-06-01")
dir.pop[,1]<-c("Ziuganov", "Mironov", "Zhirinovskii", 
               "Putin", "Mandela", "Merkel", "Nazarbaev",
               "Lukashenko", "Stalin", "Brezhnev",
               "Gorbachev", "Yeltsin", "Mikhalkov", "Grudinin",
               "Sobchak", "Kudrin","Castro","Navalny")

dir.pop[c(1:7,9:16),7]<-p.jun21
dir.pop[c(1:6,8:10,12,17),3]<-p.mar15
dir.pop[c(1:4,9:10,12),2]<-p.jan15
dir.pop[4,4]<-p.nov20
dir.pop[c(1:7,9:10,12:15,17:18),6]<-p.mar21
dir.pop[c(1:4,9:10,12:15,18),5]<-p.feb21
dir.pop[c(4,9:10,12,17),8]<-p.jun22

####Write Table for Appendix D.1
write.csv(dir.pop, file="dirpop.csv")

###Melt data and create graphic
dpop<-melt(dir.pop)
dpop[,2]<-as.Date(as.character(dpop[,2]))
names(dpop)[2]<-"Date"
dpop<-na.omit(dpop)

dpop$Figure<-as.factor(dpop$Figure)
dpop$Label<-NA
dpop$Label[which(dpop$Figure=="Putin" & dpop$Date=="2015-01-01")]<-"Putin"
dpop$Label[which(dpop$Figure=="Brezhnev" & dpop$Date=="2015-03-01")]<-"Brezhnev"
dpop$Label[which(dpop$Figure=="Castro" & dpop$Date=="2015-03-01")]<-"Castro"
dpop$Label[which(dpop$Figure=="Navalny" & dpop$Date=="2021-02-01")]<-"Navalny"
dpop$Label[which(dpop$Figure=="Grudinin" & dpop$Date=="2021-02-01")]<-"Grudinin"

dpop$Politician<-"Other"
dpop$Politician[which(dpop$Figure=="Putin")]<-"Putin"
dpop$Politician[which(dpop$Figure=="Navalny")]<-"Navalny"

dpop$Alpha<-"Other"
dpop$Alpha[which(dpop$Figure=="Putin" | dpop$Figure=="Navalny")]<-"Relevant"

ggplot(dpop,aes(Date,value,label=Label, color=Figure, alpha=Alpha)) + 
  geom_point(aes(shape=Politician), show.legend=F, size=3) + 
  geom_line(aes(linetype=Politician), show.legend=F, size=1) + scale_color_grey(start=0,end=0) + 
  geom_text(aes(color=Figure), angle=45, 
            show.legend=F,size=5, nudge_y=.05) + ylab("Proportion of Russian population \nsupporting different political figures") + 
  scale_alpha_discrete(range = c(0.25, 1)) +
  theme_bw() + coord_cartesian(ylim=c(0,1))


#################################
#####Main tables and Appendix A##
#################################

library(list)

#################
#####Table 1, A.1
#################

###2015 Jan
####Contemporary list
set.seed(123)
f21<-data.frame(cbind(jan15$pl,jan15$TreatmentP))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(putin ~ 1, data= jan15, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =  -0.057
##95% confidence interval: (-0.1643, 0.0504)

####Historical list
set.seed(123)
f21<-data.frame(cbind(jan15$hl,jan15$TreatmentH))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
predict(dim.r, direct.glm=df)
##Est. =  -0.0709
##95% confidence interval: (-0.1702, 0.0283)

###2015 March
####Contemporary list
set.seed(123)
f21<-data.frame(cbind(mar15$pl,mar15$TreatmentP))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(putin ~ 1, data= mar15, family=binomial)
predict(dim.r, direct.glm=df)
##Prediction: Difference (list - direct) 
##Est. =  -0.0841
##95% confidence interval: (-0.1927, 0.0246)

####Historical list
set.seed(123)
f21<-data.frame(cbind(mar15$hl,mar15$TreatmentH))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
predict(dim.r, direct.glm=df)
##Est. =  -0.092
##95% confidence interval: (-0.1845, 5e-04)

####Castro list (main results in text)
set.seed(123)
f21<-data.frame(cbind(mar15$il,mar15$TreatmentI))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(castro ~ 1, data= mar15, family=binomial)
predict(dim.r, direct.glm=df)
##Prediction: Difference (list - direct) 
##Est. =  -0.0876
##95% confidence interval: (-0.1916, 0.0163))

##############################
##Tables 2 (Putin 2020-1), A.1
##############################

####2020 November
####Contemporary list
set.seed(123)
f21<-data.frame(cbind(nov20$pl,nov20$TreatmentP))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(putin ~ 1, data= nov20, family=binomial)
predict(dim.r, direct.glm=df)
##Prediction: Difference (list - direct) 
##Est. =  -0.0883
##95% confidence interval: (-0.191, 0.0143)

####Historical list
set.seed(123)
f21<-data.frame(cbind(nov20$hl,nov20$TreatmentH))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
predict(dim.r, direct.glm=df)
##Est. =  -0.1372
##95% confidence interval: (-0.2289, -0.0454)

####February 2021
####Historical list
set.seed(123)
f21<-data.frame(cbind(feb21$hl,feb21$TreatmentH))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(putin ~ 1, data= feb21, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =  -0.2313
##95% confidence interval: (-0.3239, -0.1386)

####March 2021
####Historical list
set.seed(123)
f21<-data.frame(cbind(mar21$hl,mar21$TreatmentH))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(putin ~ 1, data= mar21, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =  -0.1874
##95% confidence interval: (-0.2789, -0.0959)

#### Contemporary list
set.seed(123)
f21<-data.frame(cbind(mar21$pl,mar21$TreatmentP))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
predict(dim.r, direct.glm=df)
##Est. =  -0.2351
##95% confidence interval: (-0.3387, -0.1315)

####June 2021
#### Contemporary list
set.seed(123)
f21<-data.frame(cbind(jun21$pl,jun21$TreatmentP))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(putin ~ 1, data= jun21, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =  -0.2283
##95% confidence interval: (-0.328, -0.1286)

####Intl list
set.seed(123)
f21<-data.frame(cbind(jun21$il,jun21$TreatmentI))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
predict(dim.r, direct.glm=df)
##Est. =  -0.2134
##95% confidence interval: (-0.3076, -0.1193)

###########################################################
####Tables 3 (Placebo 2021), A.3 (Placebo) and A.4 (Castro)
###########################################################

#####March
####Castro list
set.seed(123)
f21<-data.frame(cbind(mar21$il,mar21$TreatmentI))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(castro ~ 1, data= mar21, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =  -0.216
##95% confidence interval: (-0.3151, -0.1168)

#######June
####brezhnev list
set.seed(123)
f21<-data.frame(cbind(jun21$hl,jun21$TreatmentH))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(brezhnev ~ 1, data= jun21, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =   -0.2206
##95% confidence interval: (-0.3101, -0.1312)

#### Grudinin list
set.seed(123)
f21<-data.frame(cbind(jun21$fl,jun21$TreatmentF))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(grudinin ~ 1, data= jun21, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =  -0.1229
##95% confidence interval: (-0.2027, -0.043)

################################
####Tables 4 (Navalny 2021), A.2
################################

###########February
#### Navalny C list
set.seed(123)
f21<-data.frame(cbind(feb21$pl,feb21$TreatmentP))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(navalny ~ 1, data= feb21, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =  0.0127
##95% confidence interval: (-0.0822, 0.1077)

#### Navalny F list
set.seed(123)
f21<-data.frame(cbind(feb21$fl,feb21$TreatmentF))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
predict(dim.r, direct.glm=df)
##Est. =  -0.0462
##95% confidence interval: (-0.122, 0.0295)

#########March
#### Navalny F list
set.seed(123)
f21<-data.frame(cbind(mar21$fl,mar21$TreatmentF))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(navalny ~ 1, data= mar21, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =  0.0073
##95% confidence interval: (-0.0683, 0.0829)

########################################################################
###################Tables 5 (2022 results), A.1 (Putin) and A.4 (Castro)
########################################################################

####Putin historical list
set.seed(123)
f21<-data.frame(cbind(jun22$hl,jun22$TreatmentH))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(putin ~ 1, data= jun22, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =  -0.2081
##95% confidence interval: (-0.301, -0.1152)

####Putin statement list
set.seed(123)
f21<-data.frame(cbind(jun22$pi,jun22$TreatmentIP))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
predict(dim.r, direct.glm=df)
##Est. =  -0.2934
##95% confidence interval: (-0.3834, -0.2034)

####Castro intl list
set.seed(123)
f21<-data.frame(cbind(jun22$il,jun22$TreatmentI))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
df<- glm(castro ~ 1, data= jun22, family=binomial)
predict(dim.r, direct.glm=df)
##Est. =  -0.1412
##95% confidence interval: (-0.2303, -0.052)

####Castro statement list
set.seed(123)
f21<-data.frame(cbind(jun22$ci,jun22$TreatmentIC))
names(f21)<-c("hl","tH")
summary(dim.r <- ictreg(as.numeric(hl) ~ 1, data = f21, treat = "tH", J=3, method = "lm"))
predict(dim.r, direct.glm=df)
##Est. =  -0.3089
##95% confidence interval: (-0.4051, -0.2127)

################################Appendix B.1
############################Blair/Imai tests
##########Note: create objects w NAs removed

pna<-jan15$pl[which(!is.na(jan15$pl))]
tna<-jan15$TreatmentP[which(!is.na(jan15$pl))]
ict.test(pna,tna,3)

pna<-jan15$hl[which(!is.na(jan15$hl))]
tna<-jan15$TreatmentH[which(!is.na(jan15$hl))]
ict.test(pna,tna,3)

pna<-mar15$pl[which(!is.na(mar15$pl))]
tna<-mar15$TreatmentP[which(!is.na(mar15$pl))]
ict.test(pna,tna,3)

pna<-mar15$hl[which(!is.na(mar15$hl))]
tna<-mar15$TreatmentH[which(!is.na(mar15$hl))]
ict.test(pna,tna,3)

pna<-mar15$il[which(!is.na(mar15$il))]
tna<-mar15$TreatmentI[which(!is.na(mar15$il))]
ict.test(pna,tna,3)

pna<-nov20$hl[which(!is.na(nov20$hl))]
tna<-nov20$TreatmentH[which(!is.na(nov20$hl))]
ict.test(pna,tna,3)

pna<-nov20$pl[which(!is.na(nov20$pl))]
tna<-nov20$TreatmentP[which(!is.na(nov20$pl))]
ict.test(pna,tna,3)

pna<-feb21$pl[which(!is.na(feb21$pl))]
tna<-feb21$TreatmentP[which(!is.na(feb21$pl))]
ict.test(pna,tna,3)

pna<-feb21$hl[which(!is.na(feb21$hl))]
tna<-feb21$TreatmentH[which(!is.na(feb21$hl))]
ict.test(pna,tna,3)

pna<-feb21$fl[which(!is.na(feb21$fl))]
tna<-feb21$TreatmentF[which(!is.na(feb21$fl))]
ict.test(pna,tna,3)

pna<-mar21$pl[which(!is.na(mar21$pl))]
tna<-mar21$TreatmentP[which(!is.na(mar21$pl))]
ict.test(pna,tna,3)

pna<-mar21$hl[which(!is.na(mar21$hl))]
tna<-mar21$TreatmentH[which(!is.na(mar21$hl))]
ict.test(pna,tna,3)

pna<-mar21$il[which(!is.na(mar21$il))]
tna<-mar21$TreatmentI[which(!is.na(mar21$il))]
ict.test(pna,tna,3)

pna<-mar21$fl[which(!is.na(mar21$fl))]
tna<-mar21$TreatmentF[which(!is.na(mar21$fl))]
ict.test(pna,tna,3)

pna<-jun21$pl[which(!is.na(jun21$pl))]
tna<-jun21$TreatmentP[which(!is.na(jun21$pl))]
ict.test(pna,tna,3)

pna<-jun21$hl[which(!is.na(jun21$hl))]
tna<-jun21$TreatmentH[which(!is.na(jun21$hl))]
ict.test(pna,tna,3)

pna<-jun21$il[which(!is.na(jun21$il))]
tna<-jun21$TreatmentI[which(!is.na(jun21$il))]
ict.test(pna,tna,3)

pna<-jun21$fl[which(!is.na(jun21$fl))]
tna<-jun21$TreatmentF[which(!is.na(jun21$fl))]
ict.test(pna,tna,3)

pna<-jun22$hl[which(!is.na(jun22$hl))]
tna<-jun22$TreatmentH[which(!is.na(jun22$hl))]
ict.test(pna,tna,3)

pna<-jun22$il[which(!is.na(jun22$il))]
tna<-jun22$TreatmentI[which(!is.na(jun22$il))]
ict.test(pna,tna,3)

pna<-jun22$pi[which(!is.na(jun22$pi))]
tna<-jun22$TreatmentIP[which(!is.na(jun22$pi))]
ict.test(pna,tna,3)

pna<-jun22$ci[which(!is.na(jun22$ci))]
tna<-jun22$TreatmentIC[which(!is.na(jun22$ci))]
ict.test(pna,tna,3)

###########################Appendix B.2
##########################Balance tests

jun22$age<-jun22$age - mean(jun22$age)
jun22$age<-jun22$age/max(abs(jun22$age))
j22.p<-lm(TreatmentH ~ male + age + highed, data=jun22)
j22.c<-lm(TreatmentI ~ male + age + highed, data=jun22)


jun21$age<-jun21$age-mean(jun21$age)
jun21$age<-jun21$age/max(abs(jun21$age))
j21.p<-lm(TreatmentP ~ male + age + highed, data=jun21)
j21.b<-lm(TreatmentH ~ male + age + highed, data=jun21)
j21.g<-lm(TreatmentF ~ male + age + highed, data=jun21)


mar21$age<-mar21$age-mean(mar21$age)
mar21$age<-mar21$age/max(abs(mar21$age))
m21.p<-lm(TreatmentP ~ male + age + highed, data=mar21)
m21.n<-lm(TreatmentF ~ male + age + highed, data=mar21)
m21.c<-lm(TreatmentI ~ male + age + highed, data=mar21)


feb21$age<-feb21$age-mean(feb21$age)
feb21$age<-feb21$age/max(abs(feb21$age))
f21.p<-lm(TreatmentH ~ male + age + highed, data=feb21)
f21.n<-lm(TreatmentP ~ male + age + highed, data=feb21)


nov20$age<-nov20$age-mean(nov20$age)
nov20$age<-nov20$age/max(abs(nov20$age))
p.h<-lm(TreatmentH ~ male + age + highed, data=nov20)
p.c<-lm(TreatmentP ~ male + age + highed, data=nov20)


mar15$age<-mar15$age-mean(mar15$age)
mar15$age<-mar15$age/max(abs(mar15$age))
m.p<-lm(TreatmentP ~ male + age + highed, data=mar15)
m.c<-lm(TreatmentI ~ male + age + highed, data=mar15)


jan15$age<-jan15$age-mean(jan15$age)
jan15$age<-jan15$age/max(abs(jan15$age))
j.p<-lm(TreatmentP ~ male + age + highed, data=jan15)

library(stargazer)
stargazer(j.p,m.p,p.c,p.h,f21.p,m21.p,j21.p,j22.p,digits=2)
stargazer(m.c,f21.n,m21.n,m21.c,j21.b,j21.g,j22.c,digits=2)

###########################Appendix B.3
#################Ceiling/Floor analyses

###Putin historical
mean(jan15$hl[jan15$TreatmentH==F]==0,na.rm=T)
mean(jan15$hl[jan15$TreatmentH==F]==3,na.rm=T)

mean(mar15$hl[mar15$TreatmentH==F]==0,na.rm=T)
mean(mar15$hl[mar15$TreatmentH==F]==3,na.rm=T)

mean(nov20$hl[nov20$TreatmentH==F]==0,na.rm=T)
mean(nov20$hl[nov20$TreatmentH==F]==3,na.rm=T)

mean(feb21$hl[feb21$TreatmentH==F]==0,na.rm=T)
mean(feb21$hl[feb21$TreatmentH==F]==3,na.rm=T)

mean(mar21$hl[mar21$TreatmentH==F]==0,na.rm=T)
mean(mar21$hl[mar21$TreatmentH==F]==3,na.rm=T)

mean(jun22$hl[jun22$TreatmentH==F]==0,na.rm=T)
mean(jun22$hl[jun22$TreatmentH==F]==3,na.rm=T)

###Putin historical
mean(jan15$pl[jan15$TreatmentP==F]==0,na.rm=T)
mean(jan15$pl[jan15$TreatmentP==F]==3,na.rm=T)

mean(mar15$pl[mar15$TreatmentP==F]==0,na.rm=T)
mean(mar15$pl[mar15$TreatmentP==F]==3,na.rm=T)

mean(nov20$pl[nov20$TreatmentP==F]==0,na.rm=T)
mean(nov20$pl[nov20$TreatmentP==F]==3,na.rm=T)

mean(mar21$pl[mar21$TreatmentP==F]==0,na.rm=T)
mean(mar21$pl[mar21$TreatmentP==F]==3,na.rm=T)

mean(jun21$pl[jun21$TreatmentP==F]==0,na.rm=T)
mean(jun21$pl[jun21$TreatmentP==F]==3,na.rm=T)

###Putin international
mean(jun21$il[jun21$TreatmentI==F]==0,na.rm=T)
mean(jun21$il[jun21$TreatmentI==F]==3,na.rm=T)

###Putin statement
mean(jun22$pi[jun22$TreatmentIP==F]==0,na.rm=T)
mean(jun22$pi[jun22$TreatmentIP==F]==3,na.rm=T)

### Castro
mean(mar15$il[mar15$TreatmentI==F]==0,na.rm=T)
mean(mar15$il[mar15$TreatmentI==F]==3,na.rm=T)

mean(mar21$il[mar21$TreatmentI==F]==0,na.rm=T)
mean(mar21$il[mar21$TreatmentI==F]==3,na.rm=T)

mean(jun22$il[jun22$TreatmentI==F]==0,na.rm=T)
mean(jun22$il[jun22$TreatmentI==F]==3,na.rm=T)

###Statement
mean(jun22$ci[jun22$TreatmentIC==F]==0,na.rm=T)
mean(jun22$ci[jun22$TreatmentIC==F]==3,na.rm=T)

###Navalny
mean(feb21$pl[feb21$TreatmentP==F]==0,na.rm=T)
mean(feb21$pl[feb21$TreatmentP==F]==3,na.rm=T)

mean(feb21$fl[feb21$TreatmentF==F]==0,na.rm=T)
mean(feb21$fl[feb21$TreatmentF==F]==3,na.rm=T)

mean(mar21$fl[mar21$TreatmentF==F]==0,na.rm=T)
mean(mar21$fl[mar21$TreatmentF==F]==3,na.rm=T)

###Brezhnev
mean(jun21$hl[jun21$TreatmentH==F]==0,na.rm=T)
mean(jun21$hl[jun21$TreatmentH==F]==3,na.rm=T)

###Grudinin
mean(jun21$fl[jun21$TreatmentF==F]==0,na.rm=T)
mean(jun21$fl[jun21$TreatmentF==F]==3,na.rm=T)

###########################Appendix B.4
###############Direct vs list responses

#######Contemporary
ggplot(jan15,aes(politician,pl,color=as.factor(TreatmentP), fill=as.factor(TreatmentP))) + 
  geom_smooth(show.legend=F) + scale_color_grey()  + scale_fill_grey()  + coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Contemporary list, January 2015 (Putin)", x="Sum of direct control items", y="List items")

ggplot(mar15,aes(politician,pl,color=as.factor(TreatmentP), fill=as.factor(TreatmentP))) + 
  geom_smooth(show.legend=F) + scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Contemporary list, March 2015 (Putin)", x="Sum of direct control items", y="List items")

ggplot(feb21,aes(politician,pl,color=as.factor(TreatmentP), fill=as.factor(TreatmentP))) + 
  geom_smooth(show.legend=F) + scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Contemporary list, February 2021 (Navalny)", x="Sum of direct control items", y="List items")

ggplot(mar21,aes(politician,pl,color=as.factor(TreatmentP), fill=as.factor(TreatmentP))) + 
  geom_smooth(show.legend=F) + scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Contemporary list, March 2021 (Putin)", x="Sum of direct control items", y="List items")

ggplot(jun21,aes(politician,pl,color=as.factor(TreatmentP), fill=as.factor(TreatmentP))) + 
  geom_smooth(show.legend=F) + scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Contemporary list, June 2021 (Putin)", x="Sum of direct control items", y="List items")

#######Historical
ggplot(jan15,aes(historical,hl,color=as.factor(TreatmentH), fill=as.factor(TreatmentH))) + 
  geom_smooth(show.legend=F) + scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Historical list, January 2015 (Putin)", x="Sum of direct control items", y="List items")

ggplot(mar15,aes(historical,hl,color=as.factor(TreatmentH), fill=as.factor(TreatmentH))) + 
  geom_smooth(show.legend=F) + scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Historical list, March 2015 (Putin)", x="Sum of direct control items", y="List items")

ggplot(feb21,aes(historical,hl,color=as.factor(TreatmentH), fill=as.factor(TreatmentH))) + 
  geom_smooth(show.legend=F) + scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Historical list, February 2021 (Putin)", x="Sum of direct control items", y="List items")

ggplot(mar21,aes(historical,hl,color=as.factor(TreatmentH), fill=as.factor(TreatmentH))) + 
  geom_smooth(show.legend=F) + scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Historical list, March 2021 (Putin)", x="Sum of direct control items", y="List items")

ggplot(jun22,aes(historical,hl,color=as.factor(TreatmentH), fill=as.factor(TreatmentH))) + 
  geom_smooth(show.legend=F) + scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Historical list, June 2022 (Putin)", x="Sum of direct control items", y="List items")

###########Society list
ggplot(feb21,aes(figure,fl,color=as.factor(TreatmentF), fill=as.factor(TreatmentF))) + 
  geom_smooth(show.legend=F) +  scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Society list, February 2021 (Navalny)", x="Sum of direct control items", y="List items")

ggplot(mar21,aes(figure,fl,color=as.factor(TreatmentF), fill=as.factor(TreatmentF))) + 
  geom_smooth(show.legend=F) +  scale_color_grey()  + scale_fill_grey()  +coord_cartesian(ylim=c(0,4)) +
  theme_bw() + labs(title="Society list, March 2021 (Navalny)", x="Sum of direct control items", y="List items")

###########################Appendix B.5
########################List histograms

ggplot(mar15, aes(x=hl, color=as.factor(TreatmentH),
                fill=as.factor(TreatmentH))) + 
  geom_histogram(show.legend = F, position = "dodge") + 
  theme_bw() + xlab("January 2015") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(jan15, aes(x=hl, color=as.factor(TreatmentH),
                fill=as.factor(TreatmentH))) + 
  geom_histogram(show.legend = F, position = "dodge") + 
  theme_bw() + xlab("March 2015") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(nov20, aes(x=hl, color=as.factor(TreatmentH),
                  fill=as.factor(TreatmentH))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("November 2020") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(feb21, aes(x=hl, color=as.factor(TreatmentH),
                  fill=as.factor(TreatmentH))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("February 2021") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(mar21, aes(x=hl, color=as.factor(TreatmentH),
                  fill=as.factor(TreatmentH))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("March 2021") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(jun22, aes(x=hl, color=as.factor(TreatmentH),
                  fill=as.factor(TreatmentH))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("June 2022") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

####
ggplot(mar15, aes(x=pl, color=as.factor(TreatmentP),
                fill=as.factor(TreatmentP))) + 
  geom_histogram(show.legend = F, position = "dodge") + 
  theme_bw() + xlab("January 2015 (Putin)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(jan15, aes(x=pl, color=as.factor(TreatmentP),
                fill=as.factor(TreatmentP))) + 
  geom_histogram(show.legend = F, position = "dodge") + 
  theme_bw() + xlab("March 2015 (Putin)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(nov20, aes(x=pl, color=as.factor(TreatmentP),
                  fill=as.factor(TreatmentP))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("November 2020 (Putin)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(feb21, aes(x=pl, color=as.factor(TreatmentP),
                  fill=as.factor(TreatmentP))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("February 2021 (Navalny)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(mar21, aes(x=pl, color=as.factor(TreatmentP),
                  fill=as.factor(TreatmentP))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("March 2021 (Putin)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(jun21, aes(x=pl, color=as.factor(TreatmentP),
                  fill=as.factor(TreatmentP))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("June 2021 (Putin)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

###
ggplot(feb21, aes(x=fl, color=as.factor(TreatmentF),
                  fill=as.factor(TreatmentF))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("February 2021") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(mar21, aes(x=fl, color=as.factor(TreatmentF),
                  fill=as.factor(TreatmentF))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("March 2021") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()


####
ggplot(mar21, aes(x=il, color=as.factor(TreatmentI),
                  fill=as.factor(TreatmentI))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("March 2021 (Castro)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(jun21, aes(x=il, color=as.factor(TreatmentI),
                  fill=as.factor(TreatmentI))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("June 2021 (Putin)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(jun22, aes(x=il, color=as.factor(TreatmentI),
                  fill=as.factor(TreatmentI))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("June 2022 (Castro)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

####
ggplot(jun22, aes(x=pi, color=as.factor(TreatmentIP),
                  fill=as.factor(TreatmentIP))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("Putin") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(jun22, aes(x=ci, color=as.factor(TreatmentIC),
                  fill=as.factor(TreatmentIC))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("Castro") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(jun21, aes(x=pl, color=as.factor(TreatmentP),
                  fill=as.factor(TreatmentP))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("June 2021 (Putin)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(mar15, aes(x=il, color=as.factor(TreatmentI),
                fill=as.factor(TreatmentI))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("March 2015 (Castro)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(jun21, aes(x=hl, color=as.factor(TreatmentH),
                  fill=as.factor(TreatmentH))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("June 2021 (Brezhnev)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

ggplot(jun21, aes(x=fl, color=as.factor(TreatmentF),
                  fill=as.factor(TreatmentF))) + 
  geom_histogram(show.legend = F, position = "dodge")  + 
  theme_bw() + xlab("June 2021 (Grudinin)") + ylab(NULL) + 
  scale_fill_viridis_d() + scale_color_viridis_d()

##############################Appendix B.6
########################Blair/Imai Floor Effects

##############Table B.8
######Historical
f21<-data.frame(cbind(jan15$hl,jan15$TreatmentH))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
##Est. =  0.7093 (s.e. = 0.0077)

f21<-data.frame(cbind(mar15$hl,mar15$TreatmentH))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
##Est. =  0.5967 (s.e. = 0.0093) 
###Error, does not provide values on some versions of R

###Contemporary
f21<-data.frame(cbind(jan15$pl,jan15$TreatmentP))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
##Est. =  0.6659 (s.e. = 0.0083)
###Error

f21<-data.frame(cbind(mar15$pl,mar15$TreatmentP))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
##Est. =  0.6837 (s.e. = 0.0081)
###Error

#######################Table B.9
###Historical
f21<-data.frame(cbind(nov20$hl,nov20$TreatmentH))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
###Est. =  0.4865 (s.e. = 0.0104)

f21<-data.frame(cbind(feb21$hl,feb21$TreatmentH))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.3941 (s.e. = 0.0099)

f21<-data.frame(cbind(mar21$hl,mar21$TreatmentH))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.4469 (s.e. = 0.0106)
###Error

####Contemporary
f21<-data.frame(cbind(nov20$pl,nov20$TreatmentP))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
###Est. =  0.4987 (s.e. = 0.0097)

f21<-data.frame(cbind(mar21$pl,mar21$TreatmentP))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.4024 (s.e. = 0.0091)

f21<-data.frame(cbind(jun21$pl,jun21$TreatmentP))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.4376 (s.e. = 0.0095)

####International
f21<-data.frame(cbind(jun21$il,jun21$TreatmentI))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.4442 (s.e. = 0.0104)

####Table B.10
###Castro 2015
f21<-data.frame(cbind(mar15$il,mar15$TreatmentI))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
##Est. =  0.525 (s.e. = 0.0094)

####Castro 2021
f21<-data.frame(cbind(mar21$il,mar21$TreatmentI))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.3745 (s.e. = 0.0092)

####Brezhnev
f21<-data.frame(cbind(jun21$hl,jun21$TreatmentH))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.4048 (s.e. = 0.0103)

####Grudinin
f21<-data.frame(cbind(jun21$fl,jun21$TreatmentF))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.2239 (s.e. = 0.0085)
###vs 18% in list

################Table B.11
####Contemporary
f21<-data.frame(cbind(feb21$pl,feb21$TreatmentP))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.2752 (s.e. = 0.0083) 

####Figure
f21<-data.frame(cbind(feb21$fl,feb21$TreatmentF))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.1751 (s.e. = 0.0074)

f21<-data.frame(cbind(mar21$fl,mar21$TreatmentF))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.2262 (s.e. = 0.009)


######Table B.12
####Putin Historical
f21<-data.frame(cbind(jun22$hl,jun22$TreatmentH))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.5967 (s.e. = 0.0093)
###vs 63% in list

####Castro International
f21<-data.frame(cbind(jun22$il,jun22$TreatmentI))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.4095 (s.e. = 0.011)

####Castro Statement
f21<-data.frame(cbind(jun22$ci,jun22$TreatmentIC))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.2899 (s.e. = 0.0076)

####Putin Statement
f21<-data.frame(cbind(jun22$pi,jun22$TreatmentIP))
names(f21)<-c("hl","tH")
dim.r <- ictreg(hl ~ 1, data = f21, treat = "tH", J=3, method="ml",floor=T)
dim.p<-predict.ictreg(dim.r,se.fit=T)
####Est. =  0.5558 (s.e. = 0.0105)


##################Zeroes (discussed in Appendix D.6 txt)

#############################2015
##############Contemporary
mean(jan15$pl[jan15$TreatmentP==0]==0,na.rm=T) - mean(jan15$pl[jan15$TreatmentP==1]==0,na.rm=T)
###.26

mean(mar15$pl[mar15$TreatmentP==0]==0,na.rm=T) - mean(mar15$pl[mar15$TreatmentP==1]==0,na.rm=T)
###.30

##################Historical
mean(jan15$hl[jan15$TreatmentH==0]==0,na.rm=T)-mean(jan15$hl[jan15$TreatmentH==1]==0,na.rm=T)
###.20

mean(mar15$hl[mar15$TreatmentH==0]==0,na.rm=T)- mean(mar15$hl[mar15$TreatmentH==1]==0,na.rm=T)
###.30


######################################2020/1

###########Contemporary
mean(nov20$pl[nov20$TreatmentP==0]==0,na.rm=T)-mean(nov20$pl[nov20$TreatmentP==1]==0,na.rm=T)
###.14

mean(mar21$pl[mar21$TreatmentP==0]==0,na.rm=T) - mean(mar21$pl[mar21$TreatmentP==1]==0,na.rm=T) 
###.14

mean(jun21$pl[jun21$TreatmentP==0]==0,na.rm=T)-mean(jun21$pl[jun21$TreatmentP==1]==0,na.rm=T)
###.16

#########Historical
mean(nov20$hl[nov20$TreatmentH==0]==0,na.rm=T)-mean(nov20$hl[nov20$TreatmentH==1]==0,na.rm=T)
###.12

mean(feb21$hl[feb21$TreatmentH==0]==0,na.rm=T)-mean(feb21$hl[feb21$TreatmentH==1]==0,na.rm=T)
###.14

mean(mar21$hl[mar21$TreatmentH==0]==0,na.rm=T) - mean(mar21$hl[mar21$TreatmentH==1]==0,na.rm=T)
###.13

mean(jun22$hl[jun22$TreatmentH==0]==0,na.rm=T)- mean(jun22$hl[jun22$TreatmentH==1]==0,na.rm=T)
###.10



###########################Appendix C
#########################Double lists

####Jan 15 x2 list
A<-cbind(jan15$hl[jan15$TreatmentP==T],jan15$pl[jan15$TreatmentP==T]) ##Create matrix for treatment
B<-cbind(jan15$hl[jan15$TreatmentP==F],jan15$pl[jan15$TreatmentP==F]) ##create matix for control

##lw remove nas
A<-A[!is.na(A[,1]), ]
A<-A[!is.na(A[,2]), ]
B<-B[!is.na(B[,1]), ]
B<-B[!is.na(B[,2]), ]

###CIs from Glynn method
mean.jan15<- (mean(B[,1])- mean(A[,1])+mean(A[,2])- mean(B[,2]))/2
se.jan15<- sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))

mean.jan15
mean.jan15-qnorm(.975)*se.jan15
mean.jan15+qnorm(.975)*se.jan15

mean.jan15d<- mean(jan15$putin,na.rm=T)
se.jan15d <- sd(jan15$putin,na.rm=T)/sqrt(sum(!is.na(jan15$putin))) 

set.seed(123)
putin.pred.jan15 <- rnorm(10000,mean.jan15,se.jan15)
putin.pred.jan15d <- rnorm(10000,mean.jan15d,se.jan15d)

mean(putin.pred.jan15-putin.pred.jan15d)
mean(putin.pred.jan15-putin.pred.jan15d) - qnorm(.975)*sd(putin.pred.jan15-putin.pred.jan15d)
mean(putin.pred.jan15-putin.pred.jan15d) + qnorm(.975)*sd(putin.pred.jan15-putin.pred.jan15d)


######################
####March 15 x2 list
######################
A<-cbind(mar15$hl[mar15$TreatmentP==T],mar15$pl[mar15$TreatmentP==T]) ##Create matrix for treatment
B<-cbind(mar15$hl[mar15$TreatmentP==F],mar15$pl[mar15$TreatmentP==F]) ##create matix for control

##lw remove nas
A<-A[!is.na(A[,1]), ]
A<-A[!is.na(A[,2]), ]
B<-B[!is.na(B[,1]), ]
B<-B[!is.na(B[,2]), ]

###CIs from Glynn method
mean.mar15<- (mean(B[,1])- mean(A[,1])+mean(A[,2])- mean(B[,2]))/2
se.mar15<- sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))

mean.mar15
mean.mar15-qnorm(.975)*se.mar15
mean.mar15+qnorm(.975)*se.mar15


mean.mar15d<- mean(mar15$putin,na.rm=T)
se.mar15d <- sd(mar15$putin,na.rm=T)/sqrt(sum(!is.na(mar15$putin))) 

set.seed(123)
putin.pred.mar15 <- rnorm(10000,mean.mar15,se.mar15)
putin.pred.mar15d <- rnorm(10000,mean.mar15d,se.mar15d)

mean(putin.pred.mar15-putin.pred.mar15d)
mean(putin.pred.mar15-putin.pred.mar15d) - qnorm(.975)*sd(putin.pred.mar15-putin.pred.mar15d)
mean(putin.pred.mar15-putin.pred.mar15d) + qnorm(.975)*sd(putin.pred.mar15-putin.pred.mar15d)

###################################
####Feb 21 x2 list (NB: Navalny!!!)
###################################
A<-cbind(feb21$fl[feb21$TreatmentP==T],feb21$pl[feb21$TreatmentP==T]) ##Create matrix for treatment
B<-cbind(feb21$fl[feb21$TreatmentP==F],feb21$pl[feb21$TreatmentP==F]) ##create matix for control

##lw remove nas
A<-A[!is.na(A[,1]), ]
A<-A[!is.na(A[,2]), ]
B<-B[!is.na(B[,1]), ]
B<-B[!is.na(B[,2]), ]

###CIs from Glynn method
mean.feb21<- (mean(B[,1])- mean(A[,1])+mean(A[,2])- mean(B[,2]))/2
se.feb21<- sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))

mean.feb21
mean.feb21-qnorm(.975)*se.feb21
mean.feb21+qnorm(.975)*se.feb21


mean.feb21d<- mean(feb21$navalny,na.rm=T)
se.feb21d <- sd(feb21$navalny,na.rm=T)/sqrt(sum(!is.na(feb21$navalny))) 

set.seed(123)
putin.pred.feb21 <- rnorm(10000,mean.feb21,se.feb21)
putin.pred.feb21d <- rnorm(10000,mean.feb21d,se.feb21d)

mean(putin.pred.feb21-putin.pred.feb21d)
mean(putin.pred.feb21-putin.pred.feb21d) - qnorm(.975)*sd(putin.pred.feb21-putin.pred.feb21d)
mean(putin.pred.feb21-putin.pred.feb21d) + qnorm(.975)*sd(putin.pred.feb21-putin.pred.feb21d)


####################
####March 21 x2 list
####################

A<-cbind(mar21$hl[mar21$TreatmentP==T],mar21$pl[mar21$TreatmentP==T]) ##Create matrix for treatment
B<-cbind(mar21$hl[mar21$TreatmentP==F],mar21$pl[mar21$TreatmentP==F]) ##create matix for control

##lw remove nas
A<-A[!is.na(A[,1]), ]
A<-A[!is.na(A[,2]), ]
B<-B[!is.na(B[,1]), ]
B<-B[!is.na(B[,2]), ]

###CIs from Glynn method
mean.mar21<- (mean(B[,1])- mean(A[,1])+mean(A[,2])- mean(B[,2]))/2
se.mar21<- sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))

mean.mar21
mean.mar21-qnorm(.975)*se.mar21
mean.mar21+qnorm(.975)*se.mar21


mean.mar21d<- mean(mar21$putin,na.rm=T)
se.mar21d <- sd(mar21$putin,na.rm=T)/sqrt(sum(!is.na(mar21$putin))) 

set.seed(123)
putin.pred.mar21 <- rnorm(10000,mean.mar21,se.mar21)
putin.pred.mar21d <- rnorm(10000,mean.mar21d,se.mar21d)

mean(putin.pred.mar21-putin.pred.mar21d)
mean(putin.pred.mar21-putin.pred.mar21d) - qnorm(.975)*sd(putin.pred.mar21-putin.pred.mar21d)
mean(putin.pred.mar21-putin.pred.mar21d) + qnorm(.975)*sd(putin.pred.mar21-putin.pred.mar21d)

###################
####June 21 x2 list
###################
A<-cbind(jun21$il[jun21$TreatmentP==T],jun21$pl[jun21$TreatmentP==T]) ##Create matrix for treatment
B<-cbind(jun21$il[jun21$TreatmentP==F],jun21$pl[jun21$TreatmentP==F]) ##create matix for control

##lw remove nas
A<-A[!is.na(A[,1]), ]
A<-A[!is.na(A[,2]), ]
B<-B[!is.na(B[,1]), ]
B<-B[!is.na(B[,2]), ]

###CIs from Glynn method
mean.jun21<- (mean(B[,1])- mean(A[,1])+mean(A[,2])- mean(B[,2]))/2
se.jun21<- sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))

mean.jun21
mean.jun21-qnorm(.975)*se.jun21
mean.jun21+qnorm(.975)*se.jun21

mean.jun21d<- mean(jun21$putin,na.rm=T)
se.jun21d <- sd(jun21$putin,na.rm=T)/sqrt(sum(!is.na(jun21$putin))) 

set.seed(123)
putin.pred.jun21 <- rnorm(10000,mean.jun21,se.jun21)
putin.pred.jun21d <- rnorm(10000,mean.jun21d,se.jun21d)

mean(putin.pred.jun21-putin.pred.jun21d)
mean(putin.pred.jun21-putin.pred.jun21d) - qnorm(.975)*sd(putin.pred.jun21-putin.pred.jun21d)
mean(putin.pred.jun21-putin.pred.jun21d) + qnorm(.975)*sd(putin.pred.jun21-putin.pred.jun21d)

###################
####June 22 x2 list
###################

A<-cbind(jun22$pi[jun22$TreatmentH==T],jun22$hl[jun22$TreatmentH==T]) ##Create matrix for treatment
B<-cbind(jun22$pi[jun22$TreatmentH==F],jun22$hl[jun22$TreatmentH==F]) ##create matix for control

##lw remove nas
A<-A[!is.na(A[,1]), ]
A<-A[!is.na(A[,2]), ]
B<-B[!is.na(B[,1]), ]
B<-B[!is.na(B[,2]), ]

###CIs from Glynn method
mean.jun22<- (mean(B[,1])- mean(A[,1])+mean(A[,2])- mean(B[,2]))/2
se.jun22<- sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))

mean.jun22
mean.jun22-qnorm(.975)*se.jun22
mean.jun22+qnorm(.975)*se.jun22

mean.jun22d<- mean(jun22$putin,na.rm=T)
se.jun22d <- sd(jun22$putin,na.rm=T)/sqrt(sum(!is.na(jun22$putin))) 

set.seed(123)
putin.pred.jun22 <- rnorm(10000,mean.jun22,se.jun22)
putin.pred.jun22d <- rnorm(10000,mean.jun22d,se.jun22d)

mean(putin.pred.jun22-putin.pred.jun22d)
mean(putin.pred.jun22-putin.pred.jun22d) - qnorm(.975)*sd(putin.pred.jun22-putin.pred.jun22d)
mean(putin.pred.jun22-putin.pred.jun22d) + qnorm(.975)*sd(putin.pred.jun22-putin.pred.jun22d)

############################
####June 22 x2 list (Castro)

A<-cbind(jun22$ci[jun22$TreatmentI==T],jun22$il[jun22$TreatmentI==T]) ##Create matrix for treatment
B<-cbind(jun22$ci[jun22$TreatmentI==F],jun22$il[jun22$TreatmentI==F]) ##create matix for control

##lw remove nas
A<-A[!is.na(A[,1]), ]
A<-A[!is.na(A[,2]), ]
B<-B[!is.na(B[,1]), ]
B<-B[!is.na(B[,2]), ]

###CIs from Glynn method
mean.jun22c<- (mean(B[,1])- mean(A[,1])+mean(A[,2])- mean(B[,2]))/2
se.jun22c<- sqrt((1/(2*(nrow(A) + nrow(B))))*(var(A[,1]) + var(A[,2]) - 2*cov(A[,2],A[,1]))+(1/(2*(nrow(A) + nrow(B))))*(var(B[,1]) + var(B[,2]) - 2*cov(B[,2],B[,1])))

mean.jun22c
mean.jun22c-qnorm(.975)*se.jun22c
mean.jun22c+qnorm(.975)*se.jun22c

mean.jun22dc<- mean(jun22$castro,na.rm=T)
se.jun22dc <- sd(jun22$castro,na.rm=T)/sqrt(sum(!is.na(jun22$castro))) 

set.seed(123)
putin.pred.jun22c <- rnorm(10000,mean.jun22c,se.jun22c)
putin.pred.jun22dc <- rnorm(10000,mean.jun22dc,se.jun22dc)

mean(putin.pred.jun22c-putin.pred.jun22dc)
mean(putin.pred.jun22c-putin.pred.jun22dc) - qnorm(.975)*sd(putin.pred.jun22c-putin.pred.jun22dc)
mean(putin.pred.jun22c-putin.pred.jun22dc) + qnorm(.975)*sd(putin.pred.jun22c-putin.pred.jun22dc)

########################################
############Comparison of March 15 to 21
########################################

mean(putin.pred.mar15 - putin.pred.mar15d - putin.pred.mar21 + putin.pred.mar21d)
mean(putin.pred.mar15 - putin.pred.mar15d - putin.pred.mar21 + putin.pred.mar21d) + qnorm(.975)*sd(putin.pred.mar15 - putin.pred.mar15d - putin.pred.mar21 + putin.pred.mar21d)
mean(putin.pred.mar15 - putin.pred.mar15d - putin.pred.mar21 + putin.pred.mar21d) - qnorm(.975)*sd(putin.pred.mar15 - putin.pred.mar15d - putin.pred.mar21 + putin.pred.mar21d)

###############Appendix Table D.2

p.jun21<-apply(jun21[,c("ziuganov", "mironov", "zhirinovskii", "putin", "mandela",
                        "merkel", "nazarbaev", "stalin", "brezhnev", 
                        "gorbachev", "yeltsin", "mikhalkov", "grudinin",
                        "sobchak", "kudrin")], 2,
               function(x)mean(is.na(x)))
p.mar15<-apply(mar15[,c("ziuganov", "mironov", "zhirinovskii", "putin", 
                      "mandela", "merkel", "lukashenko", "stalin", 
                      "brezhnev", "yeltsin", "castro")], 2,  
               function(x)mean(is.na(x)))

p.jan15<-apply(mar15[,c("ziuganov", "mironov", "zhirinovskii", "putin",
                      "stalin", "brezhnev", "yeltsin")], 2, 
               function(x)mean(is.na(x)))

p.nov20<-mean(is.na(nov20[,c("putin")]))
p.mar21<-apply(mar21[,c("ziuganov", "mironov", "zhirinovskii", "putin", "mandela",
                        "merkel", "nazarbaev", "stalin", "brezhnev", 
                        "yeltsin", "mikhalkov", "grudinin",
                        "sobchak", "castro","navalny")], 2, 
               function(x)mean(is.na(x)))

p.feb21<-apply(feb21[,c("ziuganov", "mironov", "zhirinovskii", "putin",
                        "stalin", "brezhnev", 
                        "yeltsin", "mikhalkov", "grudinin",
                        "sobchak", "navalny")], 2, 
               function(x)mean(is.na(x)))

p.jun22<-apply(jun22[,c("putin",
                        "stalin", "brezhnev", "yeltsin", "castro")], 2, 
               function(x)mean(is.na(x)))

dir.pop<-data.frame(matrix(ncol=8,nrow=18))
colnames(dir.pop)<-c("Figure","2015-01-01","2015-03-01",
                     "2020-11-01", "2021-02-01", "2021-03-01",
                     "2021-06-01","2022-06-01")
dir.pop[,1]<-c("Ziuganov", "Mironov", "Zhirinovskii", 
               "Putin", "Mandela", "Merkel", "Nazarbaev",
               "Lukashenko", "stalin", "Brezhnev",
               "Gorbachev", "Yeltsin", "Mikhalkov", "Grudinin",
               "Sobchak", "Kudrin","Castro","Navalny")

dir.pop[c(1:7,9:16),7]<-p.jun21
dir.pop[c(1:6,8:10,12,17),3]<-p.mar15
dir.pop[c(1:4,9:10,12),2]<-p.jan15
dir.pop[4,4]<-p.nov20
dir.pop[c(1:7,9:10,12:15,17:18),6]<-p.mar21
dir.pop[c(1:4,9:10,12:15,18),5]<-p.feb21
dir.pop[c(4,9:10,12,17),8]<-p.jun22

write.csv(dir.pop, file="dirna.csv")